Mathematica for the lab - worksheet 3

Organising code and data - a trapped particle

Please don't send people copies of this worksheet. Instead, direct them to
http://j-star.org/mathematica_course.html
 where the latest version will always be available.
 © Jony Hudson, Imperial College London, 2006-2011

We're near the half-way point of the workshop, so perhaps now is an appropriate time to pause and look back on what we've covered so far.

In the first worksheet we looked at some of the features of the Mathematica language. In particular we concentrated on the so-called "higher-order" aspects of the language - functions of functions. Crucially, these functions specify what they do, not how they do it. This distinction is more important than it might seem. To give a concrete example, I could (in principle) break into your computer one evening and replace the Map function with one that dispatches individual function applications to an enormous pool of servers distributed around the world (in fact, Mathematica has a built in function ParallelMap which can do more or less that). Your code would run unchanged even though it's doing its business in a completely different way. The only change you'd notice would be that for large computations it would run much faster (given some assumptions about bandwidth!). Google uses exactly this technique to generate their web indices (or at least they did until recently) and do many of their other large computations. The beauty of it is that the same code that runs on a single workstation also runs on a (10^5) node compute farm. All because we stated what to do, not how to do it. If you're interested in Google's implementation, take a look at this paper about the MapReduce framework.

The other thing we looked at, albeit briefly, in the first worksheet was the idea of Mathematica having a uniform representation for everything. We haven't really made that much use of this yet. As we think about how to structure our programs and data more carefully this will start to become more significant.

In the second worksheet we looked at fitting and analysing data. Throughout this we focused on the role of the function as the primary mechanism for organising our code. We saw that by abstracting commonly used code into functions we could effectively extend Mathematica upwards towards our problem, eventually making the problem easy to solve. This is by no means the end of code organisation. Many of you will have encountered more sophisticated techniques for code organisation - object-oriented programming, aspect-oriented programming, <insert more buzz-phrases here>. These are all possible in Mathematica although they kind of turn out to be not as useful as you might think. There'll be a little more on this later in this sheet.

Whilst we've thought a bit about organising our code, we have yet to pay much attention to how we organise data. So far we've just used simple lists. In this worksheet we'll build a simulation of ions in a Paul trap. This will be complicated enough that it'll force us to think about how data can be structured and how code can be built around these structures. We'll see that in Mathematica, unlike in many other languages, the techniques for structuring code and structuring data have significant overlap. This is, of course, a consequence of the fact that code and data are deep-down made of the same stuff.

I should note here: this workbook gets hard towards the end, and might take quite a long time to solve. I decided it was better to present a fairly realistic example than to simplify. So don't worry too much if you get stuck or can't find the time to finish it. Hopefully the ideas will sink in, and you can look back at it sometime later if you fancy. In any case be sure to skim-read to the end as there are some juicy tidbits down there!

The Paul trap

Hopefully everybody knows roughly how a Paul trap works. The über-executive summary goes thusly: Laplace's equation tells us that in free-space a static electric field can never confine a charged particle in all directions (because one of those second derivatives must have the opposite sign to the other two to make the sum zero, so in at least one dimension the particle is not trapped). However, this is not the case for time-varying electric fields.

We'll consider a two-dimensional Paul trap which consists of  an oscillating electric field with quadrupolar symmetry. The potential, which looks like a saddle, traps the ion in one direction and forces it out in the other. We'll sinusoidally oscillate the field so that the trapped and untrapped dimensions are continually swapping. Do this right and, perhaps amazingly, the ion stays in the trap. Below I plot the potential (click the play button to start the animation).

trapPotential=Sin[(Subscript[ω, rf]) t](1/2)m ((Subscript[ω, t])^2)((x^2)-(y^2))

Animate[
Plot3D[trapPotential/.{(Subscript[ω, t])->0.1,(Subscript[ω, rf])->1,m->1,t->tt},{x,-10,10},{y,-10,10},PlotRange->{Automatic,Automatic,{-0.5,0.5}},Mesh->False],
{tt,0,2π},AnimationRate->0.6,AnimationRunning->False,SaveDefinitions->True
]

Note that I've defined my potential symbolically and then used a set of rules to substitute in "physical" values. [Note also the little trick I played with t and tt. It works around a scoping problem with Animate.]

The physics of the Paul trap is very similar to the physics of the alternating-gradient molecule decelerators that we build in the CCM lab.

Solving the equations of motion

Before we start writing any sophisticated simulation code we should make sure we can solve the equations of motion for a single ion in the trap. I'll do this as an example because there are a few little details which are a pain to get right.

The equation of motion is, of course, just Newton's law. All we need is the force on the particle as a function of time and position. To get the force we differentiate the above potential.

calculateForces[trapPotential_]:=D[trapPotential,#]&/@{x,y}/.{x->x[t],y->y[t]}
forces=calculateForces[trapPotential]

Notice a few things. First of all I'm differentiating with respect to x and y by mapping over the list {x,y}. This is no different to our previous uses of Map, but it might feel a little different because x and y don't have any values. If it's confusing remember that we're just feeding arguments to a function - nothing more. After I've built the list of the two derivatives I apply a transformation with /. that turns x into x[t] and similarly for y. This is to get the expression into a form that Mathematica's differential equation solver will be happy with (it likes to know explicitly that x and y are functions of t). This is one of the tricky bits that I mentioned above. Every time I solve a differential equation in Mathematica I still have to mess about to get this form right. Finally, I've written a function to do the calculation as we may use it again later for different trap potentials.

Now we can build the equation of motion.

makeEquationOfMotion[forces_]:=Thread[m #''[t]+γ #'[t]&/@{x,y}==forces]
eqnOfMotion=makeEquationOfMotion[forces]

Note that I've inserted a damping term proportional to the velocity - you might imagine that this coarsely models a cooling laser, or the effect of buffer gas. Also notice that I've used a shorthand for D, an apostrophe. This works because x and y are defined as functions of t only, so there's no ambiguity as to what I'm differentiating with respect to. The double equals sign which means equality, versus the single equals which means assignment, is essential. Finally note the use of Thread to force the equality signs inside the lists. An example may help to make the use of Thread clear.

{a,b}=={c,d}

Thread[{a,b}=={c,d}]

Now we're in a position to solve the equations of motion for some particular initial conditions. The Mathematica function to numerically solve differential equations is NDSolve. I'm going to pass in the initial conditions as a list of equations. Note carefully that we're solving for x[t], not for x. I find this confusing.

sol=NDSolve[
{eqnOfMotion/.{m->1,(Subscript[ω, rf])->10,(Subscript[ω, t])->2,γ->0.01},
x[0]==2,
y[0]==2,
x'[0]==0,
y'[0]==0},
{x[t],y[t]},{t,0,500},MaxSteps->∞]

NDSolve returns interpolating functions which are Mathematica's way of representing numerical solutions as functions. Internally they work by patching an approximation to the solution together from bits of polynomial (FullForm one if you want to see its guts - you might want to try a smaller one than the ones generated above). All that matters now is that if we feed them a value for their variables (in the above, t) they'll return a number. I'll plot the ion trajectory with ParametricPlot (note that I'm setting PlotPoints to higher than the default - otherwise, the plot loses a lot of detail).

ParametricPlot[{x[t],y[t]}/.sol,{t,0,500},AspectRatio->1,PlotPoints->100]

Groovy. We can see the secular motion (the big ellipses) the micromotion (the little wiggles) and the effect of the damping. That's most of ion trap physics right there!

Exercise

Solve for some different initial conditions. Try different masses and rf frequencies. Have a play! Make sure you understand the way the equations of motion were constructed. Try solving a different differential equation if you're not sure (how about the forced harmonic oscillator?).

Solutions

I've already done enough work. Over to you!

More ions, more physics, more everything

Right, so we've established that we can solve the basic problem - integrating the equations of motion in the trap. But as it stands it's a bit inconvenient and slightly useless: we have to code the initial conditions in as equations, we can only solve for one ion and we're limited to just solving for the ion rolling round the trap. We should build our code into something more solid to solve these problems.

In order to do this the first step we need to take is to think about what we actually want from the completed code. This is an often overlooked step - thinking! I’m going to spell out in detail how I’d think about this problem. It’s a bit long, but hopefully it will be useful, so bear with me. Of course, sometimes it's not possible to know what you want before you start, but usually it pays to think a bit about what sort of things you might ultimately want your code to do. Often a few simple changes/extensions at this early stage will make it much, much simpler to extend and maintain the code later on. So what would we like our code to do ? Here are a few ideas:

⋄ Simulate the trajectories of a cloud of (non-interacting) ions.
⋄ Support more than one species of ion at a time.
⋄ Allow us to track the evolution of thermodynamic quantities as the ions move in the trap (by thermodynamic I mean properties that belong to the cloud, not the individual particles).
⋄ [Advanced] Allow us to study simple interaction/collision models.

That's a lot of stuff. Where do we start?

Identifying concepts and relationships

A useful first step in designing code is to identify the distinct concepts in the thing you're modelling and try and describe the relationships between them. We can then try and capture these concepts and relationships in our data and code structures. [Note that a great deal of the hoo-ha about modern object-orientated programming as practised by people who call themselves architects is about formalising (and often providing language support for) this idea of capturing the concepts. Said architects often make a big deal out of this stuff but it's actually kind of obvious and as old as programming itself. What's new is that they've formalised it to a degree where they can demand giant salaries for doing it and then get minions to actually write the code!]

I can see the following things and concepts in the Paul trap problem:
	⋄ A particle, which has physical properies like mass, species etc.
	⋄ The position of a particle in phase space at a point in time.
	⋄ The trajectory of a particle.
	⋄ A cloud of particles.
	⋄ Thermodynamic properties of a cloud.
	⋄ The potential.
	⋄ The equations of motion.

And the relationships between these concepts? The position of a particle in phase space is associated with a particular particle. A particle can have many of these positions, each at a different time. The set of all these positions is a trajectory. Given a particle, a point in phase space, a potential, the electrode positions and equations of motion the subsequent trajectory of the particle is defined. A cloud of particles is made up of particles. The thermodynamic properties of a cloud are associated with a cloud. The equations of motion are just the equations of motion. The potential is just the potential.

The key point, that we'll explore below, is that things generally correspond to data structures, and relationships to functions.

Building the structures

Now we can try and distil these concepts and relationships into data structures and code. But before we can really do that we'll have to think just a little more about data structures.

It's becoming obvious that we're going to need a data structure to represent a particle and maybe one to represent a trajectory. How do we define these data structures? One way, the method we've implicitly used up until now, is to structure our data as lists with everything else being done by convention. By this I mean, a list {{2,3},{4,5},{6,7}} can be thought of, by convention, as a list of x, y pairs. We respect that convention when we write any function that acts on the list. This works fine for moderately small programs, but it rapidly becomes unwieldy for larger applications. For a start, it's easy to accidentally break one's own convention. And what if the data structure needs to be extended (like, in our ion trap example, we realise that the particle's internal state needs to be tracked for instance)? With simple lists inserting an extra element somewhere in the middle can often cause all of the code to break. That doesn't seem like a very maintainable or scalable approach (oooh, buzzzz!).

The real problem is that we're mixing up two different things: the abstract idea of the data structure, and the way that we actually implement the data structure. What do I mean by this? The abstract idea of a data structure could be thought of as the set of conditions that a thing has to satisfy to qualify as that kind of data structure. For our x, y pairs that's simply that we should be able to get the x- and a y-values for the ith pair. We could think about defining our data structure as anything that can return something sensible for the following functions, where I'm deliberately not even trying to write down the right-hand-side:

getX[data_,i_]:=…
getY[data_,i_]:=…

We also need a constructor function which can make a data structure.

makeXYs[listOfXYPairs_]:=…

For the implementation of this abstract data structure given above, the nested list, I can easily complete these functions.

getX[data_,i_]:=data[[i]][[1]]
getY[data_,i_]:=data[[i]][[2]]
makeXYs[listOfXYPairs_]:=listOfXYPairs

data=makeXYs[{{2,3},{4,5},{6,7}}];
getX[data,2]
getY[data,3]

It really helps to think this way, keeping distinct the ideas of the contract that defines the abstract data structure, and the particular code that implements the structure in a particular way. It may seem simple, but we've done something very important here. We've trapped the implementation of the data structure inside these three functions.

We can think of other implementations for a list of x, y  pairs. How about if I store the data as a flat list, with the convention that the data alternate x then y ?

getX[data_,i_]:=data[[2*i-1]]
getY[data_,i_]:=data[[2*i]]
makeXYs[listOfXYPairs_]:=Flatten[listOfXYPairs]

data=makeXYs[{{2,3},{4,5},{6,7}}];
getX[data,2]
getY[data,3]

In an abstract sense this is the same data structure, even though it's implemented in a different way. Notice that the code that uses the data structures is the same in both cases. Making this kind of abstraction explicit in your code is key to writing good programs. We've decoupled the rest of the program from the data-structure, making it much easier to make changes.

The trick is to try and avoid referring to the implementation of the data structure whenever possible - divide the code into the stuff that absolutely has to grub around inside the implementation of the structure and the stuff that simply needs to know that they're x, y pairs (and, as such, can always use the getX and getY functions).

To underline the utility of this idea, imagine, for some reason, I need to switch to a different implentation of x, y pairs (maybe getting them from a database because I've got so many pairs that they can't fit into my machine's memory all at once). If my code only ever gets x and y values using getX and getY then I can just rewrite the definitions of these functions and leave the rest of the code untouched. This idea is often known as "data encapsulation" or "building abstraction barriers". Once again, we're winning by specifying what to do, not how to do it.

Mathematica offers a pretty much limitless number of ways to create these data abstractions and the corresponding implementations. We'll stay focussed on the simple ones for now, and we'll come back and consider some more sophisticated structures at the end.

Mutable data

We haven't provided any way to modify the elements of the above data structure. This turns out to be a bit of a tricky point in Mathematica which we'll now try to address. [This section assumes in parts that you have fairly detailed knowledge of another programming language like C or Java - don't worry if you don't, as you probably won't have noticed that the problem I'm trying to explain is a problem!]

In most mainstream languages if we want to modify data structures passed to a function we need to use references (or pointers, which are essentially the same). By passing the data structure in "by reference" we are able to modify it and keep the modifications. The problem is that it's not very natural to pass around "references" to objects in Mathematica. [It can be done using Hold, but it's not pretty.] The easiest solution is to use a modify-reassign pattern, so our modifications will look like this:

dat=setSomeParameter[dat,newValue];

We have a function that modifies the data structure and returns the new modified data structure. We catch the new data structure and assign it back to the original symbol. This is just the same idea that you use to modify primitive types, which cannot be referenced, in C or Java:

int i = 5;
i = (3*i)+5;

The function that let's us change a part of an expression is ReplacePart. Here it is, with its bizzare syntax, used to modify an element of the data structure.

getX[data_,i_]:=data[[2*i-1]]
getY[data_,i_]:=data[[2*i]]
makeXYs[listOfXYPairs_]:=Flatten[listOfXYPairs]
setX[data_,i_,newX_]:=ReplacePart[data,(2*i-1)->newX]

And here it is in action:

data=makeXYs[{{2,3},{4,5},{6,7}}];
getX[data,3]
data=setX[data,3,9];
getX[data,3]

The particle structure

OK, with all that preparation out of the way, let's get started building the data structures. First up, I think, should be the particle. We said above that a particle has physical properties, a trajectory and a position at every time (that can be derived from the trajectory). This suggests the following sort of abstract data structure to me. We can easily add more physical properties if we like.

makeParticle[mass_,species_,trajectory_]:=…
getMass[particle_]:=…
getSpecies[particle_]:=…
getTrajectory[particle_]:=…
particle = setTrajectory[particle_,newTrajectory_]:=…

Notice that I've also included a function to change a particle's trajectory, but not the other properties. By choosing which set- and get-functions to offer I can control which of the particle's properties are mutable. I haven't included a position. My feeling is that the trajectory already contains all of the information, and we should derive positions when needed - we'll write a function to do this later.

To get you started, here’s an implementation of the above data structure using a simple list.

makeParticle[mass_,species_,trajectory_]:={mass,species,trajectory}
getMass[particle_]:=particle[[1]]
getSpecies[particle_]:=particle[[2]]
getTrajectory[particle_]:=particle[[3]]
setTrajectory[particle_,newTrajectory_]:=ReplacePart[particle,3->newTrajectory]

And here are a few examples of using it.

exampleParticle=makeParticle[5,"Ca",598];
getMass[exampleParticle]
getSpecies[exampleParticle]
getTrajectory[exampleParticle]
exampleParticle=setTrajectory[exampleParticle,thisWouldBeATrajectoryIfWeKnewHowToDefineItYet];
getTrajectory[exampleParticle]

Note that I haven’t done anything particularly clever or advanced here. People often get hung up on this bit, and you should try to avoid that! All we’re doing is giving names to the very simple functions that extract data from our lists. But by doing it we’re making our future life much easier, because we’re separating the implementation of the data structures from their logical definition.

The trajectory structure

As for the trajectory, the set of interpolating functions returned by NDSolve seems like a fine way to represent the trajectory over a certain period of time. I'm going to add some entries in the trajectory structure that say over which time period it's valid. This is convenient, but a little redundant as they're already hidden in the interpolating function. Here's what I think the structure should look like (path is the interpolating function):

makeTrajectory[tStart_,tEnd_,path_]:=…
getStartTime[trajectory_]:=…
getEndTime[trajectory_]:=…
getPath[trajectory_]:=…

Exercise

Implement this data structure. Make an example trajectory to test your functions. Just use any old set of interpolating functions (say one of the ones from the solutions above) to test the functions for now.

Solution

makeTrajectory[tStart_,tEnd_,path_]:={tStart,tEnd,path}
getStartTime[trajectory_]:=trajectory[[1]]
getEndTime[trajectory_]:=trajectory[[2]]
getPath[trajectory_]:=trajectory[[3]]

exampleTrajectory=makeTrajectory[0,500,sol[[1]]];
getStartTime[exampleTrajectory]
getPath[exampleTrajectory]

Exercise

Make an example particle with an example trajectory. Use the interpolating function generated above (in my solution example) as the trajectory's path.
Hint: note that NDSolve returns a list of solutions. You want the first (and only) interpolating function.

Solution

exampleParticle=makeParticle[40,"Ca",makeTrajectory[0,500,sol[[1]]]];

The phase space point structure

Next up is the point in phase space that a particle is occupying at a particular time. The obvious data structure for the point is:

makePhaseSpacePoint[x_,y_,vx_,vy_,t_]:=…
getX[point_]:=…
getY[point_]:=…
getVX[point_]:=…
getVY[point_]:=…
getT[point_]:=…

Exercise

I'm sure you're getting the idea by now. Implement the above functions.
Note: you might be thinking that is all getting rather tedious. You'd be right. Maybe we could write some code to do it for us ...

Solution

makePhaseSpacePoint[x_,y_,vx_,vy_,t_]:={x,y,vx,vy,t}
getX[point_]:=point[[1]]
getY[point_]:=point[[2]]
getVX[point_]:=point[[3]]
getVY[point_]:=point[[4]]
getT[point_]:=point[[5]]

examplePointInPhaseSpace=makePhaseSpacePoint[10,5,0,0,1000];

Capturing relationships

We've built our data structures. Now let's capture some of the relationships between them as functions.

Exerise

Write a function that extracts a point in phase space, using the functions defined above, from a trajectory, given a time. Define your function as pointFromTrajectory[trajectory_,time_]:=…. Test your function on your example particle's trajectory.
Hint: you can differentiate interpolating functions just like anything else.
Hint 2: you might find the function's a little bit fiddly. By dealing with the fiddliness here, we don't have to worry about it later.
Hint 3: note that the critical point is that you shouldn’t access any of the data structures directly in this code. You should use the encapsulated methods defined above (i.e. getX[point], not point[[1]]).

Solution

pointFromTrajectory[trajectory_,time_]:=Module[{p,vals},
p=getPath[trajectory];
vals={p,D[p,t]}/.t->time;
makePhaseSpacePoint[x[time],y[time],x'[time],y'[time],time]/.Flatten[vals]
]

pointFromTrajectory[getTrajectory[exampleParticle],50]

Exercise

Make the trivial but convenient extension to pointFromParticle[particle_,time_]:=….

Solution

pointFromParticle[particle_,time_]:=pointFromTrajectory[getTrajectory[particle],time]

pointFromParticle[exampleParticle,50]

Solving the equation of motion

We can finally get on to the meat of the problem, which now turns out to be very easy: solving the equation of motion.

Exercise

Write a function that takes a particle, an initial point in phase space (as defined above), a set of parameters to substitute into the equation of motion, a duration to solve for (i.e. it solves the equations from the time in the starting point to the time in the starting point plus the duration to solve for) and a maximum number of steps. It should return the particle with the resulting trajectory attached.

Make your function definition look like solveInTrap[particle_,startingPoint_,parameters_,Δt_,maxSteps_]:= …. Test your function with the arguments below. Solve for t between 0 and 20. Pass the options PrecisionGoal->3,AccuracyGoal->3 to NDSolve.
Hint: note that NDSolve returns a list of solutions, which should only be one element long as the solution should be unique. That means you want to take the first part of what NDSolve returns.
Hint 2: once again, this is a function that modifies the particle, so you should catch the modified particle that it returns in an assigment.

p=makeParticle[1,"Ca",{}];
sp=makePhaseSpacePoint[10,5,0,0,0];
parameters={(Subscript[ω, rf])->10,(Subscript[ω, t])->2.1,γ->0};

Solution

solveInTrap[particle_,startingPoint_,parameters_,Δt_,maxSteps_,potential_:trapPotential]:=Module[{startTime,sol},
startTime=getT[startingPoint];
sol=NDSolve[
{
makeEquationOfMotion[calculateForces[potential]]/.m->(getMass[particle])/.parameters,
x[startTime]==getX[startingPoint],
y[startTime]==getY[startingPoint],
x'[startTime]==getVX[startingPoint],
y'[startTime]==getVY[startingPoint]
},
{x[t],y[t]},
{t,startTime,startTime+Δt},
MaxSteps->maxSteps,
PrecisionGoal->3,
AccuracyGoal->3
];
setTrajectory[particle,makeTrajectory[startTime,startTime+Δt,sol[[1]]]]
]

p=solveInTrap[p,sp,parameters,20,∞]

Time for a break ...

That might seem like quite a lot of work to achieve more or less nothing. But that's not the case. By starting to capture the ideas and relationships between the concepts underlying the program, and by abstracting away the implementation of the data structures, we'll find that subsequent coding becomes much easier. Again, it's a case of building up layer by layer until we have a language for solving the problems we want to solve. We've built the foundations, so let's make something of them.

Back again. Viewing the solutions

Before we get too much further let's write a function that plots the trajectory of a particle. Write a function to plot the trajectory between two times. And put a dot on the particle’s position at the end time. Include a parameter which lets you set the plot range, pr - you should plot from {-pr,pr} in both dimensions. Make sure your plot function does not try and plot for values of t outside the range you have solved for.

Test your function by plotting the trajectory of the particle you solved for above between t=5 and t=15. A plot range of 20 should be appropriate.

Solution

plotParticleTrajectory[particle_,startTime_,endTime_,plotRange_]:=Module[{endPoint},
plotStart=Max[startTime,getStartTime[getTrajectory[particle]]];
plotEnd=Min[endTime,getEndTime[getTrajectory[particle]]];
endPoint=pointFromParticle[particle,plotEnd];

Show[
ParametricPlot[
{x[t],y[t]}/.getPath[getTrajectory[particle]],{t,plotStart,If[plotEnd==plotStart,plotEnd+(10^(-6)),plotEnd]},AspectRatio->1,PlotPoints->200,PlotStyle->Blue
],
Graphics[{Red,PointSize[0.02],Point[{getX[endPoint],getY[endPoint]}]}],
PlotRange->{{-plotRange,plotRange},{-plotRange,plotRange}}
]
]

Note that I've taken care of the case where plotStart and plotEnd are the same.

plotParticleTrajectory[p,5,500,20]

Make an animation

For a bit of fun, make an animation of the ion going round in the trap, say from t=0 to t=20. If you're feeling fancy, density plot the time-varying potential in the background as well. Admire the magic of the Paul trap.

Solution

plotTrapPotential[parameters_,tEnd_,plotRange_,zNorm_]:=DensityPlot[
((trapPotential/.m->1/.parameters/.t->tEnd)/zNorm),
{x,-plotRange,plotRange},
{y,-plotRange,plotRange},
PlotRange->{Automatic,Automatic,All},
Mesh->False,
ColorFunction->GrayLevel,
ColorFunctionScaling->False
]

plotWithPotential[particle_,parameters_,tStart_,tEnd_,plotRange_]:=Show[
plotTrapPotential[parameters,tEnd,plotRange,600],
plotParticleTrajectory[particle,tStart,tEnd,plotRange],
PlotRange->{{-plotRange,plotRange},{-plotRange,plotRange}}
]

Animate[plotWithPotential[p,parameters,0,t,20],{t,0,20},AnimationRate->0.8,AnimationRunning->False]

Clouds of particles

Let's look at clouds of ions. For the moment we're going to completely neglect the interaction between the ions (much more like a molecule decelerator than an ion trap). I think here we probably don't need a complicted data structure: we can represent a cloud of particles as a list of particle structures.

Generate a cloud of 50 particle phase space positions centred on (10,10) with zero average velocity, give the points Gaussian position and velocity distributions of width 5 and 1 respectively at t=0. Your result should be a list of phase space positions as defined above. We'll be making more than one cloud, so make this into a function, makePhaseSpaceCloud.

noOfParticles=50;

Solution

makePhaseSpaceCloud[n_,t_,{x_,σx_},{y_,σy_},{vx_,σvx_},{vy_,σvy_}]:=
Table[makePhaseSpacePoint[RandomReal[NormalDistribution[x,σx]],RandomReal[NormalDistribution[y,σy]],RandomReal[NormalDistribution[vx,σvx]],RandomReal[NormalDistribution[vy,σvy]],t],{n}]

initialCloudPositions=makePhaseSpaceCloud[noOfParticles,0,{10,5},{10,5},{0,1},{0,1}];

Simulating the cloud's motion

Solve the equations of motion for these initial conditions, with the following list of particles. Solve using the parameters below. Solve up to t=100.

particles=Table[makeParticle[1,"Ca",{}],{noOfParticles}];
parameters={(Subscript[ω, rf])->10,(Subscript[ω, t])->2.1,γ->4 (10^(-2))};

Solution

particles=MapThread[solveInTrap[#1,#2,parameters,100,∞]&,{particles,initialCloudPositions}];

A function to display the cloud

Write a function to plot the positions of the particles in a cloud at time t. Your function should take a plot range parameter as before.

Plot the cloud at some time. Make an animated plot of the cloud's motion up until t=50. Try putting the potential in the background again.

Solution

plotCloud[cloud_,parameters_,t_,plotRange_,showTrapPotential_:True]:=
Show[
If[showTrapPotential,plotTrapPotential[parameters,t,plotRange,3000],Graphics[]],
ListPlot[
{getX[#],getY[#]}&[pointFromParticle[#,t]]&/@cloud,PlotStyle->{Red,PointSize[0.015]}
],
PlotRange->{{-plotRange,plotRange},{-plotRange,plotRange}}
]

plotCloud[particles,parameters,20,50]

Animate[plotCloud[particles,parameters,t,50],{t,0,50},AnimationRunning->False,AnimationRate->1.5]

Thermodynamic properties of the cloud

Size

Write a function to extract the cloud's size at time t. Plot the cloud size as a function of time between t=0 and t=100.

Solution

sizeOfCloud[particles_,time_]:=StandardDeviation/@Transpose[({getX[#],getY[#]}&[pointFromParticle[#,time]]&/@particles)]

sizes=sizeOfCloud[particles,#]&/@Range[0,100,0.1];

ListPlot[Transpose[sizes],Joined->True]

Temperature

Write a function to calculate the temperature of the cloud at time t. Plot the temperature as a function of t for the cloud above.
Hint: you might have to think a bit about what temperature means. In particular, does it even make sense to talk of temperature in a trap? Is the cloud in thermal equilibrium? What's the most sensible thing to derive?

Solution

temperatureOfCloud[particles_,time_]:=StandardDeviation[(Sqrt[Plus@@{((getVX[#])^2),((getVY[#])^2)}])&[pointFromParticle[#,time]]&/@particles]

temps=temperatureOfCloud[particles,#]&/@Range[0,100,0.2];

avgTemps=Mean/@Partition[temps,{10}];

ListPlot[avgTemps,Joined->True,PlotRange->All]

Multi-species clouds

Let's simulate a cloud with more than one species of ion. We will find that because of the thought we put in at the start that this is very straightforward. In particular we can assign the particles a color and track this through the simulation without any extra work.

Extend the data structure

First of all, patch the data structure to support a new physical property, color. You'll need to rewrite the constructor and add a new "get" function for the color. Test your functions.

Solution

makeParticle[mass_,species_,trajectory_,color_:Red]:={mass,species,trajectory,color}
getColor[particle_]:=particle[[4]]

pcl=makeParticle[1,"Mg",{},Blue];
getColor[pcl]

Generate a mixed cloud

Generate a cloud that has two species of ion one called "Ca" with mass 1 and another called "Mg" with mass 0.6. Generate 50 of each ion. Asssign each particle a color, Blue for Ca and Red for Mg.

Solution

numberOfEachSpecies=50;
initialCloudPositions2=makePhaseSpaceCloud[2 numberOfEachSpecies,0,{10,5},{10,5},{0,1},{0,1}];

particles2=Join[Table[makeParticle[1,"Ca",{},Blue],{numberOfEachSpecies}],Table[makeParticle[0.6,"Mg",{},Red],{numberOfEachSpecies}]];

Simulating the cloud's motion

Solve the equations of motion for these particles. Solve using the parameters below. Solve up to t=100.

parameters={(Subscript[ω, rf])->10,(Subscript[ω, t])->2.1,γ->4 (10^(-2))};

Solution

particles2=MapThread[solveInTrap[#1,#2,parameters,100,∞]&,{particles2,initialCloudPositions2}];//AbsoluteTiming

Plotting the cloud

Write a function to plot the cloud. Color the particles appropriately. Again, if you're feeling fancy show the oscillating potential in the background. Try to make your solution general for any number of colors if possible. I mean, try not to have a Red and a Blue in there. You might find this is quite tricky!

Solution

plotMultiSpeciesCloud[cloud_,parameters_,t_,plotRange_]:=
Module[{colors,subClouds,plots},
colors=Union[getColor/@cloud];
subClouds[c_]:=Select[cloud,getColor[#]==c&];
subClouds=subClouds/@colors;
getPlotPoints[subCloud_]:={getX[#],getY[#]}&[pointFromParticle[#,t]]&/@subCloud;
plots=MapThread[
ListPlot[getPlotPoints[#1],PlotStyle->{#2,PointSize[0.015]}]&,
{subClouds,colors}
];
Show@@(Join[
{plotTrapPotential[parameters,t,plotRange,3500]},plots,
{PlotRange->{{-plotRange,plotRange},{-plotRange,plotRange}}}])
]

Animate[plotMultiSpeciesCloud[particles2,parameters,t,50],{t,0,50},AnimationRunning->False,AnimationRate->1.5]

Pseudopotential

A common approximation is to ignore the micromotion - the little wiggles - and just concentrate on the secular motion in the trap. You can show that the secular motion is well described by a "pseudopotential" of the form

pseudopotential=-(1/2)m ((ωpp)^2)((x^2)+(y^2))

where the secular trapping frequency, ωpp, is given by ((Subscript[ω, t])/((Sqrt[2])(Subscript[ω, rf]))). Let's try modelling the ion behaviour in this pseudopotential approximation.

Extend the data structure

Go back and modify your solveInTrap function so that it accepts a potential as a parameter. Make sure that this doesn't break the rest of the code that you've written so far. Simulate two particles, one in the real potential and one in the pseudo potential using the parameters below. What do you see?

p=makeParticle[2,"Ca",{}];
p2=makeParticle[20,"Ca",{}];
sp=makePhaseSpacePoint[2,2,20,-30,0];
parameters={(Subscript[ω, rf])->100,(Subscript[ω, t])->10 (Sqrt[(Sqrt[2])]),γ->0};
psParameters={ωpp->1,γ->0};

Solution

p2=solveInTrap[p2,sp,parameters,20π,∞];
p=solveInTrap[p,sp,psParameters,20π,∞,pseudopotential];

Show[plotParticleTrajectory[#,0,4π,50]&/@{p,p2}]

Looks good. There's still some divergence from the initial micromotion kick that the real particle gets.

Parametric heating

Thus far we've assumed an ideal trap. How about if we include some instability?

We might imagine that our trap potential shakes due to some technical noise. Perhaps the power supplies that supply the voltage to the electrodes are imperfect. We can model this with a new (pseudo)potential:

pseudopotentialShaking=-(1/2)m ((ωpp)^2)(1+δ)((x^2)+(y^2))

Since the ion trap is harmonic, we might imagine that it's useful to look at the response of the system to sinusoidal perturbations of the trap strength.

Model a cloud of ions in the imperfect trap. Shake the trap sinusoidally at a frequency of your choosing. Take a look at what the ions are doing as a function of the shaking frequency. You might find the parameters below provide a useful starting point.

noOfParticles=100;
initialCloudPositions3=makePhaseSpaceCloud[ noOfParticles,0,{0,20},{0,20},{0,200},{0,200}];

Solution

psParameters={ωpp->10,γ->0.1,δ->0.2 Sin[40t]};
particles3=Table[makeParticle[1,"Ca",{},Blue],{noOfParticles}];

particles3=MapThread[solveInTrap[#1,#2,psParameters,20,∞,pseudopotentialShaking]&,{particles3,initialCloudPositions3}];//AbsoluteTiming

Animate[plotCloud[particles3,parameters,t,50,False],{t,0,20},AnimationRunning->False,AnimationRate->2]

Frequency dependence

Let's examine the frequency response more closely. Write a function that models a cloud of ions in the trap, letting them evolve for a certain amount of time with a certain shaking frequency, and then calculates their temperature. Look at how the temperature depends on frequency.
Hint: the interesting frequency is 2 ωpp.
Note: you might find that this is tremendously memory hungry. The reason is subtle and annoying. It can be worked around with a little effort!

Solution

tempAfterParametricHeating[noOfParticles_,psParameters_,tMax_]:=Module[{particles,particles2,initialCloudPositions},

particles=Table[makeParticle[1,"Ca",{},Blue],{noOfParticles}];
initialCloudPositions=makePhaseSpaceCloud[ noOfParticles,0,{0,20},{0,20},{0,200},{0,200}];
particles2=MapThread[solveInTrap[#1,#2,psParameters,tMax,∞,pseudopotentialShaking]&,{particles,initialCloudPositions}];
temperatureOfCloud[particles2,tMax]
]

tempAfterParametricHeating[200,{ωpp->10,γ->0.1,δ->0.1 Sin[20.1t]},20]

parametricTemps={#,tempAfterParametricHeating[50,{ωpp->10,γ->0.1,δ->0.1 Sin[#t]},20]}&/@Range[15,25,0.05];//AbsoluteTiming

ListPlot[parametricTemps,Joined->True,PlotRange->All]

Looks ok. There's a resonant response around 2 ωpp. The structure on the peak comes from not waiting long enough for the cloud to equilibrate.

Coulomb crystals (advanced - you might want to skip it for now)

So far we've completely neglected the interactions between the ions. Some fun stuff happens in interacting ion clouds. One particularly pretty thing is that if the cloud is cooled enough then the ions arrange themselves into a regular "crystal". Let's see if we can model this.

We'll need to rewrite some of the fundamentals of our simulator to make this work, so I present it as an open exercise. What needs to be changed? What can we re-use? What, in retrospect, should we have written differently above?

Some suggestions: try using a trap which has different x- and y-trapping frequencies; try and model clouds with more than one type of ion; how about giving different ion species different damping constants to simulate sympathetic cooling, where one species is cooled and cools other species in the trap through the Coulomb interaction.

Solution

Here's my solution. It's far from perfect. It's pretty far from nice for that matter. I decided to change as little as possible above, as those solutions would be confusing otherwise. But if I was going to change the above functions I think I'd get rid of the x[t] and y[t] in the trajectory data structure, as they are particularly awkward. I'd probably extend the force and equation of motion functions to take coordinate symbols as parameters. And I'd modify the earlier graphing routines to respect the ions' colors.

Still, maybe this solution will be interesting because it shows some (in some cases quite nasty) tricks that you can use to work yourself out of a hole!

squashedPseudoTrap=-(1/2)m(((ωx)^2) (x^2)+((ωy)^2) (y^2));

I'm going to build the coupled equations of motion fully symbolically.

coulombForce[{x1_,y1_},{x2_,y2_}]:=(1/((((((x1-x2))^2)+(((y1-y2))^2)))^(3/2))){x1-x2,y1-y2}

forceOnParticle[otherParticlePositions_,potential_]:=calculateForces[potential]+Plus@@(coulombForce[{x[t],y[t]},#]&/@otherParticlePositions)

forceOnParticle[{{x2,y2},{x3,y3}},squashedPseudoTrap]

makeSymbol[base_,n_]:=Symbol[base<>ToString[n]]

equationOfMotionForInteractingParticle[particleIndex_,noOfParticles_,potential_]:=Module[{otherParticleRange,makeParticleCoords,otherParticleCoords},
otherParticleRange=DeleteCases[Range[noOfParticles],particleIndex];
makeParticleCoords:={makeSymbol["x",#][t],makeSymbol["y",#][t]}&;
otherParticleCoords=makeParticleCoords/@otherParticleRange;
makeEquationOfMotion[forceOnParticle[otherParticleCoords,potential]]/.Thread[{x,y,γ,m}->{
makeSymbol["x",particleIndex],
makeSymbol["y",particleIndex],
makeSymbol["γ",particleIndex],
makeSymbol["m",particleIndex]
}]
]

equationOfMotionForInteractingParticle[1,3,squashedPseudoTrap]

makeCoupledEquationsOfMotion[noOfParticles_,potential_]:=Flatten[equationOfMotionForInteractingParticle[#,noOfParticles,potential]&/@Range[noOfParticles]]

makeCoupledEquationsOfMotion[2,squashedPseudoTrap]

Remove[makeParticle];
makeParticle[mass_,species_,trajectory_,color_:Red,γ_:0]:={mass,species,trajectory,color,γ}
getγ[particle_]:=particle[[5]]

There's some crufty symbol manipulation in here that could be avoided if made the modifications I described above.

solveInTrapInteracting[particles_,startingPoints_,parameters_,Δt_,maxSteps_,potential_]:=Module[{startTime,reps,initialConditions,solveFor,sol},
startTime=getT[startingPoints[[1]]];
reps=Flatten[{makeSymbol["m",#]->getMass[particles[[#]]],makeSymbol["γ",#]->getγ[particles[[#]]]}&/@Range[Length[particles]]];
initialConditions=Flatten[{
makeSymbol["x",#][startTime]==getX[startingPoints[[#]]],
makeSymbol["y",#][startTime]==getY[startingPoints[[#]]],
makeSymbol["x",#]'[startTime]==getVX[startingPoints[[#]]],
makeSymbol["y",#]'[startTime]==getVY[startingPoints[[#]]]
}&/@Range[Length[particles]]];
solveFor=Flatten[{makeSymbol["x",#][t],makeSymbol["y",#][t]}&/@Range[Length[particles]]];
sol=NDSolve[
Join[
makeCoupledEquationsOfMotion[Length[particles],potential]/.reps/.parameters,
initialConditions
],
solveFor,
{t,startTime,startTime+Δt},
MaxSteps->maxSteps,
PrecisionGoal->3,
AccuracyGoal->3
];
MapThread[setTrajectory[#1,makeTrajectory[startTime,startTime+Δt,{x[t]->#2[[1]][[2]],y[t]->#2[[2]][[2]]}]]&,{particles,Partition[sol[[1]],{2}]}]
]

Again, I could avoid writing this if I modified the plotting functions above.

plotCloudSimple[cloud_,parameters_,t_,plotRange_,showTrapPotential_:True]:=
Show[
Graphics[
MapThread[{getColor[#1],Disk[{getX[#2],getY[#2]},0.3]}&,{cloud,pointFromParticle[#,t]&/@cloud}]
],
PlotRange->{{-plotRange,plotRange},{-plotRange,plotRange}}
]

numberOfEachSpecies=20;
initialCloudPositions4=makePhaseSpaceCloud[ 2 numberOfEachSpecies,0,{0,5},{0,5},{0,5},{0,5}];
particles4=Join[Table[makeParticle[1,"Ca",{},Blue,0],{numberOfEachSpecies}],Table[makeParticle[0.6,"Mg",{},Red,0.5],{numberOfEachSpecies}]];

particles4=solveInTrapInteracting[particles4,initialCloudPositions4,{ωx->1,ωy->2},180,∞,squashedPseudoTrap];//AbsoluteTiming

Animate[plotCloudSimple[particles4,{},t,20,False],{t,0,180},AnimationRunning->False,AnimationRate->3]

Cool!

Other data structures

So far we've been using simple lists as the implementation of our data structures. These are very straightforward to use, but can become unweildy for complex structures. Mathematica doesn't offer a library of data structures like most programming environments do, but it's rather striaghtforward to roll your own. In this section we'll look at a couple of structures that might prove useful.

Rule-based associative arrays

Associative arrays are found in almost all programming environments (known variously as hashmaps, dictionaries, lookup tables ...). They are lists of key, value pairs. Look them up on Wikipedia if you're not familiar with the idea. The abstract set of rules that define an associative array might be (with the hopefully obvious associated semantics):

makeAssociativeArray[{keyValuePairs_}]:=…
getValue[aa_,key_]:=…
aa=setValue[aa_,key_,newValue_]:=…
aa=addEntry[aa_,key_,value_]:=…

Notice that I've included aa=… in the specification. As before this is to indicate that the function returns the modified data structure.

In Mathematica associative arrays can be very easily implemented with rules. Here's an implementation of the above function specifications.

makeAssociativeArray[keyValuePairs_List]:=Rule@@@keyValuePairs
getValue[aa_,key_]:=key/.aa
setValue[aa_,key_,newValue_]:=aa/.(key->_)->(key->newValue)
addEntry[aa_,key_,value_]:=Append[aa,key->value]

Aren't they beautifully concise? Make sure that you understand how they work. If in doubt, close this worksheet and try and write your own. In particular make sure you understand the setValue function. The underscore stands for "anything". It's the same expression you've been using when defining functions, although that might not be obvious right now!

Here is an associative array in action.

testAA=makeAssociativeArray[{{a,1},{b,2},{c,3}}]

getValue[testAA,b]

testAA=setValue[testAA,c,5];
getValue[testAA,c]

testAA=addEntry[testAA,d,"fingers"]

getValue[testAA,d]

It's a very powerful and flexible data structure.

Exercise

Implement a function removeEntry[aa_,key_] that allows removal of entries from the array. Your function should return the new array with the entry removed. Test your function by removing the entry for d in the above test array. Don't forget that you need to assign the result back to testAA if you want to keep the change.

Solution

removeEntry[aa_,key_]:=DeleteCases[aa,(key->_)]

testAA=removeEntry[testAA,d]

I don't know what to call it

Another neat way to store data by key is to make use of Mathematica's pattern matcher. We can define what looks like a function, and give it values for "special cases":

dat[1]=100;
dat[5]=1000;

We can then read out the values of this "function" to retrieve the data:

dat[5]

So, this behaves just as an associative array, as we described above. We're not restricted to numeric values for either the argument or the result.

dat2["name"]="hey, steve";
dat2["age"]=100;
dat2["lucky numbers"]={1,4,76,19};

You can take a look at what values are associated with a symbol with the ? function (a.k.a. Information).

?dat2

Conclusion

We've covered a lot in this worksheet. We've learned how to solve differential equations. We've thought about designing data structures and learned techniques for building them. This only scratches the surface of program design, although it does cover a lot of what's needed to be productive.

The key thing that we've been doing is building abstraction barriers. In our code above we have at least three identifiable levels of abstraction. At the highest level there's code that only needs to know that particles, points, trajectories exist - things like our plotting functions and solving function. Beneath that there's the code that implements the points, particles, trajectories that needs to know that lists exist. And then finally there's the lowest level that implements low-level operations on lists that's built in to Mathematica. [It doesn't stop there of course, there are many more levels before we get down to machine code!]

By layering the code like this we've made it very flexible. We could decide that we no longer want to use flat lists to implement the data structures. All we'd have to do is rewrite a handful of fairly trivial functions. In fact ...

Exercise

You've decided, for some reason or other, that you no longer like flat lists. Re-implement the ion trap data structures in a new way. Make sure that your simulation code runs unchanged by simulating a cloud and plotting the results. Feel the awesome power of abstraction!

Solution

Ok, so just to make the point that we're really completely isolated from the underlying implementation of the particle data structure by the formal interface I'm going to implement it with the most unusual data structure I can think of. Here are three functions (the weird names come from LISP) that can be used to build an access lists. Notice that the data is stored as functions - try and figure out how the functions work and what they do.

cons[x_,y_]:=Function[m,m[x,y]]
car[z_]:=z[Function[{p,q},p]]
cdr[z_]:=z[Function[{p,q},q]]

l=cons[3,cons[4,cons[5,{}]]]
car[l]
car[cdr[l]]
car[cdr[cdr[l]]]

I'll re-implement the particle using these funny lists.

makeParticle[massIn_,speciesIn_,trajectoryIn_]:=cons[massIn,cons[speciesIn,cons[trajectoryIn,{}]]]
getTrajectory[particle_]:=car[cdr[cdr[particle]]]
setTrajectory[particle_,newTrajectory_]:=particle/.(getTrajectory[particle]->newTrajectory)
getMass[particle_]:=car[particle]
getSpecies[particle_]:=car[cdr[particle]]

Quick check that it works.

p=makeParticle[5,"Ca","weeee"]
getMass[p]
getSpecies[p]
getTrajectory[p]
p=setTrajectory[p,"woooo"];
getTrajectory[p]

And the other structures.

makeTrajectory[tStartIn_,tEndIn_,pathIn_]:=cons[tStartIn,cons[tEndIn,cons[pathIn,{}]]]
getStartTime[trajectory_]:=car[trajectory]
getEndTime[trajectory_]:=car[cdr[trajectory]]
getPath[trajectory_]:=car[cdr[cdr[trajectory]]]

makePhaseSpacePoint[xIn_,yIn_,vxIn_,vyIn_,tIn_]:=cons[xIn,cons[yIn,cons[vxIn,cons[vyIn,cons[tIn,{}]]]]]
getX[point_]:=car[point]
getY[point_]:=car[cdr[point]]
getVX[point_]:=car[cdr[cdr[point]]]
getVY[point_]:=car[cdr[cdr[cdr[point]]]]
getT[point_]:=car[cdr[cdr[cdr[cdr[point]]]]]

Solve as before, with the same code!

p=makeParticle[1,"Ca",{}];
sp=makePhaseSpacePoint[10,5,0,0,0];
parameters={(Subscript[ω, rf])->10,(Subscript[ω, t])->2.1,γ->0};
p2=solveInTrap[p,sp,parameters,20,∞]
plotParticleTrajectory[p2,5,15,20]

Magic isn't it. The power of abstraction!

Pragmatism

Before we close this section on thinking about what we've done I feel a final note is in order. To be honest with you, most of the code I write doesn't use formalised, abstract data structures like we've used above. Most stuff I write in Mathematica is sufficiently small and hackish that I just use lists or lists of rules directly. There's a balance to be struck between getting everything perfect and well organised and just getting it done. It's up to you to find that balance in your own code. I almost always think about the code I'm writing in this "abstract data structure" way though, even if I don't write it like that. And certainly, for larger programs I abstract the data structures explicitly.

Postcript: types

In the last section of the worksheet we'll take a quick look at a topic related to structuring data which I think is fun and interesting.

Mathematica doesn't seem to worry much about what type of data it's dealing with. This is in stark contrast to languages like C, Java etc where you have to tell the compiler exactly what sort of data to expect for every variable and function (double, int, char etc). Often, especially for the sort of problems that one frequently solves in Mathematica, this can be a great advantage - it certainly makes coding a lot less painful.

But there are advantages to having types. Typing is often used in the "other" languages to build data structures and enforce rules about their usage, as well as placing constraints on function behaviour. At the simplest level, structs can be defined in C to group primitive data together. The object-oriented languages (Java, C++, C# etc) let you define objects which are, in some sense, a type of structured data. The key thing is that these languages let you define new types and make behaviour conditional upon those types. Can we do this in Mathematica?

Well, yes and no. Yes, because, of course, we can do anything in Mathematica, and no because it might not look quite like what we're used to. The technique is very simple, but if you think about it too much it becomes mind-numbingly complicated. The basic idea is that we can define our own types by making expressions with a symbol that denotes the type at their head. For instance I can define a type, human, as follows.

jony=human[35,84,"Brown"]

Now, you might be thinking that this just looks like a call to a function that I haven't defined yet. And indeed it does. But that's just because you're used to interpreting such expressions as code. They can just as well be a kind of data. I'm never going to define a function "human"; it's a marker for a type. People familiar with algebraic data types might call the above a constructor (this might confuse object-oriented programmers though as they tend to think of constructors as actual function calls!).

What's the point of this? For starters I can restrict functions to only accept data of type human.

hairColor[h_human]:=h[[3]]

hairColor[jony]
hairColor[{30,69,"Brown"}]
hairColor["Steve"]

So this is slightly useful, but the real value comes when we specialise the behaviour of existing functions with respect to our newly defined type. This gets a bit complicated to understand, because it requires knowledge of the way Mathematica evaluates expressions that you probably don't have yet. It's pretty simple to do though, so let's look at an example. We'll worry about how it works later.

human/:Print[h_human]:=Print["Age: "<>ToString[h[[1]]]<>", Hair color: "<>hairColor[h]", Weight: "<>ToString[h[[2]]]]

Print[jony]

I define what looks like a normal function prefixed by "human /:". This tells Mathematica to associate the definition with the symbol human and not with Print as it normally would. The Mathematica way of saying this is that the definition is an upvalue of human, rather than a downvalue of Print as it otherwise would be.

It's probably too much to get into what all of this means at the minute, but you'll notice that we've associated the behaviour of a function with a particular type of data. This is starting to look like the beginnings of object-oriented programming. [In fact, it's pretty straightforward to build a full-blown object system on top of Mathematica if you need it]. For the object-oriented programmers amongst you, note that operator overloading is possible in exactly the same way. Let's define addition for humans as making a human with the average of their ages, the sum of their weights, and blue hair.

human/:Plus[h1_human,h2_human]:=human[((h1[[1]]+h2[[1]])/2),h1[[2]]+h2[[2]],"Blue"]

monster=jony+human[41,55,"Green"]
Print[monster]

This is quite likely all a bit confusing. I include it here just to show what's possible. And also because I think it's interesting to think about the very blurred distinction between code and data.

A final note: you've been using these kinds of data structures from the beginning of the first worksheet, and probably from the very first time you used Mathematica. Both Graphics objects and complex numbers are implemented this way. The fact that it probably hadn't struck you as strange that Graphics objects look like functions that have no action is testament to how powerful a technique they are for structuring and encapsulating code and data. [Try FullForm on a graph or a complex number.]

Exercise

If you're not sick of it all yet, try rewriting the ion simulator using types as defined above!

Solution

You're on your own with this one!